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Abstract 

Kepler's Equation is solved over the entire range of elliptic motion by a 
fifth-order refinement of the solution of a cubic equation. This method is not 
iterative, and requires only four transcendental function evaluations: a 
square root, a cube root, and two trigonometric functions. The maximum 
relative error of the algorithm is less than one part in 10 18 , exceeding the 
capability of double-precision computer arithmetic. Roundoff errors in 
double-precision implementation of the algorithm are addressed, and 
procedures to avoid them are developed. 

Introduction 

Numerical solution of the two-body problem for orbital motion is heavily dependent on 
efficient solution of Kepler's Equation 


M — E — esin E 


( 1 ) 


for the eccentric anomaly E in terms of the eccentricity e and mean anomaly M [ 1 , 2]. Most 
methods involve choosing a starting formula and then improving this using an iterative refinement 
method. Many of these methods have difficulties in the critical region where eccentricity is close to 
one and mean anomaly is close to zero; some iterative methods even fail to converge in this region 
[3]. There have been several comparisons of numerical methods for solving Kepler's Equation, 
with conflicting claims for the accuracy and efficiency of the various algorithms [3-8]. 

This paper presents a new algorithm using a starting formula resulting from solution of a cubic 
equation based on a Pad<S approximation to the sine function. This starting formula has smaller 
errors than any previously considered [3-8]. A single application of a fifth-oider method is used to 
refine the starting estimate, rather than iteration of a lower-order method to satisfy a convergence 
criterion. The latter procedure would require at least two more trigonometric function evaluations, 
and perhaps many more. Odell and Gooding, among others, have emphasized the advantages of 
refinement using a fixed number of iterations [3]. As pointed out by Mikkola [7], such a method is 
not really iterative, but is actually a direct solution of Kepler's Equation to the desired accuracy. 

The method has errors that are less than the least significant bit of double-precision floating- 
point numbers over the entire range of elliptic motion, 0 < e < 1. The limit of unit eccentricity is not 
really elliptic motion, and the solution of Kepler's Equation is not useful there, but the eccentricity 
can be arbitrarily close to this limit. We present a derivation of the algorithm and display contour 
plots of its errors as functions of the eccentricity and eccentric anomaly. Numerical problems 
arising in double-precision implementation of the algorithm are also discussed and resolved. 


47 





Starting Formula 


Our method starts with a Padd approximation for sin E, depending on a parameter a : 


sin E ~ o(a,E) = 


6a-(ft-3)E 2 
6a + 3E 2 


( 2 ) 


It is assumed that E and M have been reduced by multiples of 2/r to have absolute value less than 
or equal to n. The Taylor series expansion of this approximation at E = 0 is 


£ 3 E 5 


(3) 


This expansion is exact through terms of order E 3 , which is crucial for good performance in the 
critical region with e near unity and M near zero. The series is exact througl^ terms of order E 3 for 
a = 10, and the approximation for sin E is exact at E- + 7rfor a = 3n /(n -6) =7.65. The 
precise specification of the parameter cc for our method will be considered below. 

Substitution of equation (2) into Kepler's Equation gives the cubic equation 


Defining 

and 


[3(1 -e)+ ae]E 3 - 3 ME 2 + 6a(l - e)E - 6aM = 0. 
d = 3(1 — e) + cte 

z = dE 

j2 • 


and multiplying through by d 2 gives the standard form for the cubic equation [9] 


with 


z 3 + a^z +a\z + ciQ — 0, 

«2 = —3 M , 


and 


a\ = bad(\-e), 
oq = -6otd 2 M . 


(4) 

(5) 

( 6 ) 

(7) 

(8a) 

(8b) 

(8c) 


We define the auxiliary quantities 


q = a l /3-aj/9 = 2ad(\- e)-M 2 

(9) 

= (a\a 2 - 3ag )/6 - a\ /27 = 3ad(d - 1 + e)M + M 3 

(10) 


and note that equation (7) has a unique real root if q 3 + r 2 > 0. It is not difficult to see that this 
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condition is satisfied for ail positive a , which covers the range of interest, except in the case that 
e = 1 and M = 0, in which case the root z = 0 is threefold degenerate. In any case, the desired root 
of equation (7) is given by 


where 


and 


z = s x +s 2 -a 2 l?> = s l +, 

Si 4- M 

(id 



1 


s \ = 

l^r + ijq 3 + r 2 

r 

(12a) 

1 

n 2 \ 

1 

3 


S 2 S \ 

r-yq + r J 


(12b) 


The accurate evaluation of equation (1 1) for small M requires cancellations between s\ and si. We 
avoid this numerical problem by employing a trick attributed to Karl Stumpff by Battin [2], Write 


, _ +*2)(*l 2 -^2+^) 5 3 +j 3 2r 

+ s 2 ~ 2 ~2 -_ 2 ~ 2 ~ ~2 2 - 0 3 ) 

■*1 “ S|^2 + s 2 ^1 -^52+^2 •*! - V[.V2 + si 


We next multiply the numerator and denominator by 


w = 



(14) 


and use the fact that s x s 2 = -q to simplify this expression. The absolute value of r is used in 
equation (14) to avoid cancellations of positive and negative quantities. Inserting the result into 
equation (1 1) and using equation (6) gives the first-order solution to Kepler's equation 


El 



2rw 

+ wq + q 2 



(15) 


Since r is proportional to M as the latter quantity goes to zero, this expression does not depend on 
cancellations in that limit. The solution is singular when both e - 1 and M = 0 simultaneously, but 
this is not of concern since equation (1) shows that E = 0 whenever M = 0. Thus a numerical 
solution of Kepler's Equation is unnecessary in this case. 

This solution is the starting formula for our algorithms. It requires two transcendental function 
evaluations; the square root and cube root in equation (14). The cube root actually involves two 
transcendental function evaluations if it is implemented as a logarithm and an exponential. 
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Specification of a 


The criterion for choosing a is the minimization of the relative errors in the starting formula. 
These errors are computed on a grid of 201 values of e between 0 and 1 and 25 1 values of E 
between 0 and n. The exact value of M is calculated at each grid point from equation (1), and then 
Ej is obtained from equation (15). The relative error at each grid point is then 

error = (E,-E)/E. (16) 

These errors were computed using quadruple precision floating point numbers with 1 12 bits in the 
mantissa. The contours of constant errors for a = 10 and a = 37T 2 / (/r - 6) are shown in Figures 
1 and 2, respectively. The contours are linearly spaced with an increment of 0.001 between 
contours. The errors in Figure 1 are all negative, with a minimum value of - 0.040. The errors in 
Figure 2 are all positive and significantly smaller in magnitude, having a maximum value of 0.013. 
The errors in both figures are monotonically increasing functions of the eccentricity, which 
suggests that an optimal a could be found by minimizing the errors for eccentricity equal to unity. 


-0.040 



Figure 1. Relative Errors in Starting Formula E\ for a - 10 
Equal -error contours with 0.001 linear contour spacing 
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It is not necessary that a be a constant parameter; we can choose it to be a function of e and M. 
Equation (2) would be exact if a could be chosen to satisfy 


3£ 2 (£-sin£ ) 

a ideal “ „3 • 

E -6(£-sin£) 


(17) 


This equation clearly not useful as it stands, since E is not known until Kepler's Equation has been 
solved. However, it is possible to find a usable function a(e,M ) that is a good approximation to 
equation (17). The right side of this equation is an even function of E, and so a(e,M) must be an 
even function of M. It is also useful to choose a form that takes the value a = 37T 2 / ( n 2 - 6) when 
E = M = ± n, since this will assure the continuity of the solutions over the extended range of these 
variables outside [ - k, n ]. This is not really important, since we intend to refine the result to the 
full accuracy of machine arithmetic, but Figure 2 shows that this condition is likely to lead to a 
good starting formula. 



eccentric anomaly (degrees) 


Figure 2. Relative Errors in Starting Formula E\ for a = 3 7T 2 / ( n 2 - 6) 
Equal-error contours with 0.001 linear contour spacing 
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As discussed above, it is most important to find an accurate form for a(e,M) on the boundary 
e = 1. Figure 3 shows aj deal and the straight line fit 


a(l,M) 


37T 2 +0.87r(;r-|M|) 


iF- 6 


(18) 


where M is given by equation (1) with e = 1. This formula can be extended to all eccentricities by 
noting that equations (17) and (18) are functions of E and M, respectively. We obtain the correct 
dependence of the slope of the straight line at the right end point by making use of 


dM 

dE 


I E=n 


l + e. 


(19) 


The resulting approximation for all eccentricities and mean anomalies is given by 


a(e , M) jj^Apm±A 

7T -6 


( 20 ) 
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Our starting formula is given by the solution of the cubic for this form of a. The error contours for 
this starting formula are shown in Figure 4, with linear contour spacing of 2 x 10 - 5 . This is Fifty 
times finer than the spacing of the contours in Figures 1 and 2. The values of these errors are 
between - 2.3 x 10 and 2.8 x 10 ~ 4 , which are smaller than the errors of any other starting 
formula known to the author. The starting formula does involve a fair amount of computation, but 
Mikkola [7] has also proposed a starting formula that requires solution of a cubic equation and has 
maximum errors seven times as large as those of the method proposed here. 

Refinement 

The iterative refinements are based on finding roots of polynomial approximations of 

f(E) = E-esinE- M. (21) 

The aim of this paper is to find a computational method yielding errors are smaller than the least 
significant bit of double-precision floating point numbers with 52 bits in the mantissa, which is 



Figure 4. Relative Errors in Starting Formula E\ for a Given by Equation (20) 
Equal-error contours with 2 x 10 ~ 5 linear contour spacing 
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about 10 _ 16 . A fifth-order refinement of the starting formula is expected to have relative errors 
less of (2.8 x 10 -4 ) 5 = 1.7 x 10 ~ 18 , which is adequate to achieve double-precision accuracy. 
The third-order (Halley) and higher-order corrections are given by 


HE) = -~ 


m 


s 4 (e) = -- 


f'(E)-±f(E)f"(E)/f'(Ey 

m 


and 


S 5(E) = -~ 


/'<£)+ i< 5 3 (£)/"(£) + 5*3 (£)/"'(£) 

/(£) 


/'(£) + ^ S 4 (E)f"(E) + 1 «!(£)/'"(£) + i Sl(E)f"(E) ’ 


and 


( 22 ) 

(23) 

(24) 


where the subscripts denote the order of the correction. The required partial derivatives are: 

/'(£) = l-<? cos £, (25) 

/"(£) = esin£. (26) 

/"'(£) = ecos£ = 1 - /'(£), (27) 

f""(E) = -esinE = -f"(E). (28) 


where the second forms of equations (27) and (28) are used to avoid additional trigonometric 
function evaluations. The fifth-order refined estimate is given by 


£ 5 = £l+« 5 (£,). 


(29) 


The equal-error contours of £5 are plotted in Figure 5 with logarithmic spacing of the contours, 
which is to say that the errors on adjacent contours differ by a factor of ten. The zero contour is not 
plotted because of numerical roundoff; this contour is not significant since the magnitude of the 
errors is more important than their sign. The magnitude of the relative errors is less than 10 _24 
over large areas of the figure, and the maximum error magnitude over the entire range of elliptic 
motion is 7.35 x 10 ~ 19 . This is about half the naive prediction based on the order of the correction 
employed. This method requires only two trigonometric function evaluations in addition to the 
transcendental function evaluations need to solve the cubic equation for the starting formula. 

Numerical Considerations 


The errors of our method are smaller than the least significant bit of double-precision floating 
point numbers. Special care must be taken when this method is actually implemented in double- 
precision arithmetic, however. A naive implementation yields unacceptably large errors exceeding 
5x10 ~ 14 . Plotting error contours shows that all errors with magnitudes in excess of 5 x 10 -16 
are in the region e > 0.75 and £ < 45°. This effect was discussed by Odell and Gooding [3], who 
attribute it to cancellations in the computation of /(£) and /'(£) by equations (21) and (25), 
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7.35e-19 



Figure 5. Relative Errors in Refined Eccentric Anomaly Eg 
Equal-error contours with x 10 logarithmic contour spacing 

respectively. The solution to this problem is to modify the computation of these quantities. The 
revision to the first derivative computation is straightforward; equation (25) is replaced by 

f'(E) = l — e + 2csin 2 (£/2). (30) 

The fix to equation (21) is not so simple. The method employed here is to replace equation (21) by 
the equivalent form 

f(E) = M*(e,E)-M , (31) 

where 

M*{e,E) = E-esmE (32) 

outside the range e > 0.5 and E\ < 1 radian. Inside this range, M*(e,E) is given by the Pad6 
approximant: 

M\e,E) = (l-e)E + eE 3 n J^^L , (33) 

denominator 
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with 


numerator s -1.7454287843856404 x KT 6 £ 6 + 4.1584640418181644 x 10 E 4 

. . (34a) 

- 3.0956446448551 138 x 10“ 2 £" + 1 
and 

denominator = 1.78043671 195 19884 x 1(T 8 E 8 +5. 97276 1373 1070647 x KT 6 # 6 

(34b) 

+ 1. 0652873476684 1 42 x 1 0~ 3 E 4 + 1. 1 426 1 32 1 308693 1 7 x 1 0“ 1 E 2 + 6. 

The second derivative is computed as 

f"(E) = E-M*(e,E) (35) 

over the entire range of eccentricity and eccentric anomaly. The third and fourth derivatives are 
computed as usual. No additional transcendental function evaluations are required by these fixes, 
although the Padd approximant is probably more expensive to evaluate than the sine function. 

The resulting double precision implementation of the algorithm has relative errors less than 
4 x 10 _16 over the entire range of elliptic motion. Inspection of a plot of the error contours 
revealed no systematic pattern, confirming that the errors are due solely to unavoidable machine 
arithmetic roundoff. 


Discussion 

The algorithm developed in this paper has been shown to have errors well within the inherent 
limitations of double-precision computer arithmetic, over the entire range of elliptic orbital motion. 
Among algorithms with this property, this method is at least as efficient as any proposed 
previously, requiring only four transcendental function evaluations: a square root, a cube root, and 
two trigonometric functions. The method is singular only when the eccentricity is unity and the 
mean anomaly is simultaneously zero. There are two reasons why it is never necessary to solve 
Kepler’s Equation in this case: first, unit eccentricity is not really elliptic motion, and second, the 
eccentric anomaly is known to be zero when the mean anomaly is zero, making numerical solution 
unnecessary. Special procedures to handle double-precision roundoff errors near this singular 
point have been developed and tested. The resulting algorithm is well suited for implementation in 
orbit propagation systems. 
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